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During the past six or seven years, several methods of orbit refinement 
were developed specifically for use with artificial satellites and spacecraft. 
This article describes these methods, and the classical method, in a uniform 
mathematical formalism in order to facilitate comparisons of their relative 
advantages and disadvantages for practical systems applications. However, 
such comparisons are made in this article only to the extent that motivated 
the development of the new methods. 

I. INTRODUCTION 

The accuracy to which the position of a satellite or spacecraft can be 
predicted depends upon the accuracy to which the "initial" position and 
velocity vector, or related orbital parameters, arc known. Since these 
parameters can be determined only from observational data which in- 
evitably contain observational errors, the accuracy to which they can be 
known depends upon the nature, the quantity, the accuracy, and the dis- 
tribution (in space and time) of the observational data, and the way in 
which these data are processed. The accuracy of the orbital parameters, 
and of prediction, depend also upon the accuracy to which all of the forces 
acting on the satellite or spacecraft are known and taken into account. 
Clearly, the term "accuracy" must be taken largely in a statistical sense. 

Orbit refinement is essentially data smoothing for the purpose of ac- 
curate prediction. Given the nature of the observational data, and the 
statistical properties of the observational errors, it is possible to formu- 
late a method of data smoothing and prediction which is optimum in the 
sense of giving predictions with the greatest possible accuracy. However, 
such an optimum method will, in general, be useful only as a standard of 
comparative performance for more practical methods. The reason for 
this is that it has not been difficult to find simpler and therefore more 
practical methods which are nearly as accurate as the optimum method. 
(For example, see Ref. 1.) 

It should be noted also that in the practical applications of data- 
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smoothing and prediction methods each application is usually character- 
ized by a unique set of practical constraints. Thus, it has been a frequent 
experience that a method which was judged to be the most practical for 
one application was usually not judged to be the most practical for an- 
other application. In fact, it has frequently happened that a new method, 
perhaps new only in the sense that it is a composite of parts of older 
methods, was developed for a particular system. 

The classical method of orbit refinement, the so-called "differential 
corrections" method (which is essentially the method of least squares, 
developed by K. F. Gauss in 1795) has served astronomers very well for 
over 150 years. However, this method becomes quite unwieldy for large 
quantities of observational data. Hence, the quantity of data to be proc- 
essed was frequently reduced by the supplementary use of "normal 
places" and/or "smoothing," described briefly in Section 6.6.1, pp. 141- 
142 of Ref. 2, and later in this article. Thus, up to about 1955 no need 
was felt for another method of orbit refinement, although the basis for 
the development of alternative methods was implicit in a number of 
articles published in the field of general statistical analysis, such as Refs. 

3-6. 

With the development of artificial satellites and space probes, the 
need for alternative methods of orbit refinement began to be felt in 
some quarters. The first definite proposal of an alternative method, as 
far as the author is aware, was made by P. Swcrling (Refs. 7-9). A some- 
what different method, independently developed by the author (Ref. 10) 
was used in the Telstar I experimental satellite communications system 
(Refs. 11, 12). Some difficulties experienced with this method after about 
four weeks of successful operation led A. J. Claus (Ref. 13) to develop 
another method which is slightly different from Swerling's method. In 
addition to these methods, it is worthwhile to include a method of space 
navigation described by R. II. Battin (Ref. 14) because it involves a prac- 
tical detail which, under favorable circumstances, may be profitably 
introduced into the other methods. 

The essential details of these methods will be described here in a uni- 
form mathematical formalism in order to reveal their basic similarities 
and differences, and in order to facilitate comparisons of their relative 
advantages and disadvantages for practical systems applications. 

II. CLASSICAL DIFFERENTIAL CORRECTIONS METHOD. LEAST SQUARES 

Let (p be an n-rowed vector representation of the observational (angu- 
lar) data. Every component of this vector is assumed to be labeled to 
identify it as either a declination angle or a right ascension angle, and to 
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specify the time at which it was observed. Let e be the ti-rowed vector 
representation of a set of values of the orbital elements, and let ^(e) be 
an n-rowed vector representation of the angles which would have been 
observed, assuming; that the actual orbital elements are exactly repre- 
sented by e, and assuming that observations are made with ideal ac- 
curacy. If the observational errors are independently and normally dis- 
tributed, with zero means and equal variances, the best estimate of the 
orbital elements is that value of e which minimizes the quadratic form 

where the prime stands for transposition. This quadratic form is simply 
the sum of the squares of the components of the vector difference 
Q> — <p(e). 

Let e be an initially assumed value for e, close to the true value. Then, 
to the first-order term in (e — e ), 

p(c) = <p(to) + J(e - eo), (1) 

where J is the n by G matrix symbolized by 

./ = d V (e )/deo. (2) 

Hence, to second-order terms, 

Q = [ r - ./.( e - €„)]'• [/• - J-{e - e„)], 

where /• is the /(-rowed vector residual 

r = £-?(«o). (3) 

Now, Q is a minimum with respect to e if 

J'-[r - J(e - eo)] = 0. 

Written in the form 

./'../• (e - eo) = ./'•/-, 

this corresponds to the set of (i equations commonly called "normal equa- 
tions." If ./'•./ is nonsingular, and if the value of e which satisfies this 
equation is denoted by e, then, 

e = e n + (J'-J)~ l -J'-r. (4) 

This e is then substituted for eo in ( 2 ), in (3), and in the right-hand mem- 
ber of (4), in order to obtain another e. This substitution procedure is 
iterated until e has essentially converged. The final e is the least squares 
estimate of the orbital elements. 
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In a statistical sense, this estimate is unbiased to the first order in the 
observational errors, assuming that the errors are not biased. It is only 
asymptotically unbiased to the second order, but it is a "consistent" 
estimate in the sense that the probability is unity that it will be correct 
to the second order as n — > ». (See Section 4.3 for clarification.) 

2.1 Classical Method. Weighted Least Squares 

The method described in the preceding section is quite unwieldy for 
large values of n on account of the size of the J matrix. Some relief was 
obtained by resorting to various artifices. For example, if a number of 
observations were made at sufficiently short intervals of time, a straight 
average of these observations would be taken and treated as a single 
observation called a "normal place." (More elaborate methods of deriv- 
ing normal places directly from the observational data are called 
"smoothing" by Baker and Makemson in Ref. 2.) Since normal places 
would be more accurate than single actual observations, in proportion 
to the number of actual observations which went into each of them, it 
was necessary to generalize the differential corrections method to some 
extent. 

The quadratic form to be minimized is now 

Q= [<p-<p(e)Y-W-[<p-<p(e)}, 

where W is an n by n diagonal matrix. In expanded form, it is 

Q = T.wa-[ipi - <Pi(e)f, 

where the Wu are the components of W. Thus, the quadratic form is 
simply a weighted sum of the squares of the components of the vector 
difference £ — <p(e). 

The analysis in this case will not be pursued beyond this point, since 
it is a special case of the analysis given in the next section. Suffice it to 
say that if the analysis were carried out for this special case, the results 
would be equivalent to the method used by astronomers when they deal 
with normal places, or with uncorrelated observations of different de- 
grees of accuracy. 

2.2 General Form of the Classical Method. 

If the observational data are not all of the same nature (angles, ranges, 
and range-rates), and especially if some of the errors in the data are cor- 
related, let $ be the n by n co variance matrix of the w-rowed vector <p. 
For further generality, let e be an ra-rowed vector representation of a set 
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of values of parameters which, in addition to the G orbital elements, may 
include such quantities as the frequency of a satellite-borne Doppler 
source, and instrumental biases. Then the quadratic form to be mini- 
mized is 

Q = I* -rtO]'-*"" 1 -^ -*(«)]• (5) 

Under the assumption that the components of Cp obey a joint w-dimen- 
sional normal (i.e., Gaussian) distribution with covariance matrix <f>, 
the value of e which minimizes the quadratic form (5) is the maximum 
likelihood estimate of the true value of the parameters. Under any other 
symmetrical probability distribution of the errors in the observational 
data, this value of e is simply the weighted least squares estimate of the 
parameters. 

Substituting (1) into (5) we get 

Q = [ r - J-(e - »)Y+-*'W ~ /•(• - eo)], (G) 

where ./ and r are defined by (2) and (3) except that J is now an n by 
m matrix, and e is an m-rowed vector. Now, Q is minimum with respect 
to eif 

./'•*-■[/• -./•(«- *)] = 0. 

Denoting the value of e which satisfies this equation by e, we have 

= «o + C-p, (7) 



e 



where 



C= (./'•<!> -J) , an m by m matrix, (8) 

p = J' •$~'-r, an m-rowed vector. (9) 

Equation (7) is the generalization of (4) and, as in the case of (4), it 
is to be used iteratively until I has essentially converged. 

After e has converged, C is its covariance matrix. This follows from 
the fact that (6) may be expressed in the form 

Q = (e — e)'-C _1 -(e — e) + terms independent of e. (10) 

In case the data are all of the same nature (all angles, or all ranges, or 
all range-rates), are all of the same accuracy, and the errors are not 
correlated, then, C = a 2 -(J' J )~ l . 

The availability of the covariance matrix C of the estimate e offers the 
possibility of using the classical method in the intrapass stage of the pass- 
by-pass method described in Section IV, in order to reduce the amount 
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of observational data to be processed at any one time. However, it is not 
the only way to reduce the amount of observational data to be processed 
at any one time. There are other ways which, in particular, avoid the 
computation of the n by m matrix J, where n may be of the order of 400 
for each pass. One such way is Swerling's method described in the next 
section. Another such way is cited as an example in Section 4.3. A modi- 
fied form of the second way is described briefly in Section V. 

in. swerling's method 

Let ei be a 6-rowed vector estimate of the osculating orbital elements 
at epoch h , and let (\ be its covariance matrix. Let $ be the n-rowed 
vector of new observational data more or less centered at epoch t 2 , 
where t 2 > k , and let <f> be its covariance matrix. To obtain the least 
squares estimate of the osculating orbital elements at epoch t 2 , we must 
first update (i.e., extrapolate) «i and Ci . If ei is the result of updating 
ei , the updated Ci is 

& = M-CvM', 

where M is the 6 by 6 matrix symbolized by 

M = dei/dei . 

Then, assuming that the errors in the new data are not correlated with 
the errors in the old data, the quadratic form to be minimized is 

Q = ( e - ei)'-Cr l -(e - k) + ]S> - K01'"*"" 1 '!* ~ *(01- (H) 

This is essentially the sum of the right-hand members of (5) and (10). 
Now, to the first-order term in (e — ci), 

?(e) = *(M) +./•(* - h), (12) 

where ./ is the n by 6 matrix symbolized by 

./ = dv(h)/dk. (13) 

Then, to second-order terms, 

Q= ( e - e,)'-Cr'-(e- h) + [r - J ■ (e - h)]'-*'* ' V ~ J- (c - h)], 

where 

r = ip — <p(ei). (14) 

Now, Q is minimum with respect to e if 

A~*-(i - k) - J''*"*'\r-J-(* - k)] = o. 



METHODS OF ORBIT REFINEMENT 891 

Denoting the value of e which satisfies this equation by e, we have 

i = 6, + C-p, (15) 

where 

c = [cr 1 + ./'•*-' -./]-', (i6) 

p = ./'.$-'.,•. (17) 

Finally, since the quadratic form, to second-order terms, may be ex- 
pressed in the form 

Q = (e — l)'-G~ «(e — e) + terms independent of e, (18) 

it follows that C is the covariance matrix of e. Equations (11), (15), 
(16), and (17) correspond respectively to equations (16), (25), (19), 
and (23) in Swerling's JAS paper (Ref. 8). Further, since ci and C x may 
have been computed at the preceding stage in exactly the same Avay as 
e and C at the last stage, (11) corresponds also to equation (30) in 
Swerling's JAS paper. 

With regard to the updating of ei and C\ it should be noted that esti- 
mates of "time of perigee (or nodal) passage" should be labeled with 
the serial number of the passage. Thus, even in the hypothetical case of 
a purely Keplerian orbit, unless the serial number of the passage is in- 
tended to be the same in e as it is in ei , the vector «i will differ from the 
vector ei . The component T x (time of perigee or nodal passage) of d 
will be increased to t x in e x , where 7\ is T\ plus an integral multiple of 
the period estimate 27rai 3/2 /V/z, where ai is the semimaj or axis component 
of ei . The matrix .1/ will therefore be a unity matrix except for an off- 
diagonal component dti/dai which is an integral multiple of 3ttV ' a\i '/i. 
In this connection, A. J. Claus has found advantages in using the period 
instead of the semimajor axis as a component of «i , especially when per- 
turbing forces are taken into account. 

In the classical method n must be at least equal to the number of 
orbital elements, and it must include all of the available observational 
data — old observational data (which has been processed at least once 
before) as well as new. In Swerling's method n may be less than the num- 
ber of orbital elements (possibly n = 1), and old observational data are 
represented by the 6 components of ei , the 21 distinct components of the 
symmetrical matrix Cj , and the epoch d . The chief objection to Swerl- 
ing's method, as far at least as some applications are concerned, is its 
inability to omit any part of the old observational data without reproc- 
essing the remainder of it. This objection is less serious if Swerling's 
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method is used only in the intrapass stage of the pass-by-pass method 
described in the next section. 

IV. PASS-BY-PASS METHOD 

This method consists essentially of two stages per pass. (These stages 
are not quite the same as Swerling's stages.) In the first (or intrapass) 
stage, a set of estimates of the orbital elements, and an associated co- 
variance matrix, are computed from the observational data for each pass. 
In this stage, any method of computation which yields a covariance 
matrix for the estimates of the orbital elements may be used. In the sec- 
ond (or interpass) stage, the sets of single-pass (i.e., intrapass) estimates 
of the orbital elements are combined cumulatively (and possibly selec- 
tively), by a method in which the single-pass (i.e., intrapass) covariance 
matrices play important roles. 

Let ei and C\ have the same significance as in the description of Swerl- 
ing's method, but let the new data be processed separately to obtain an 
independent vector estimate e 2 of the osculating orbital elements at 
epoch k , with covariance matrix C 2 . Then, the quadratic form to be 
minimized is 

Q= (e- eO'-Cr'-^ " 6.) + (e - e 2 ) / -C 2 " 1 -(e - «,). (19) 

This is the sum of two terms similar to the right-hand member of (10). 
Now, Q is minimum with respect to e if 

A -1 -(i- *) + cr 1 -(«- «) = o. 

Denoting the value of « which satisfies this equation by e, we have 

- e = c-(cr 1 -6 1 + c 2 - 1 - f2 ) (20) 

where 

C = (ft" 1 + Cf 1 )-'. (21) 

Finally, since (19) may be expressed in the form 

Q = ( e - €)'-C _1 -(e - e) + terms independent of e, (22) 

it follows that C is the covariance matrix of e. (See Appendix for a more 
illuminating derivation). 

Note that a "fading memory" can be introduced into the pass-by-pass 
method (or any other method in which the old data are represented by 
ei , d , and h) simply by substituting kC t for Ci , where A; is a scalar 
constant (k > 1) or 

k = exp [y(k — t\)] 
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where 7 is a constant (7 > 0). However, at the cost of storing a number 
of single-pass estimates of the orbital elements and the associated co- 
variance matrices, old single-pass estimates of the orbital elements may 
be completely omitted at any time without having to reprocess any of 
the observational data on which the more recent single-pass estimates 
are based. 

4.1 An Example 

It is illuminating to see how the method of combination described in 
the preceding section works out in a simple problem. Consider an object 
traveling along the x axis at a known acceleration a. Let 



ei = 



be the estimates of position and velocity at epochs t\ and t 2 , and let 

* - * - [i j} 

If t 2 — h = t, then 

Vxi + vn + far 2 "] 
ei_ L ft + ir J' 

d«l ["I 7"! 

[2,22 2-1 

0-. +T* 9 T*. I 
t<j v a v J 



whence 



Then, 



Note that the correlation coefficient, which is 

1 

is very close to unity for large values of t. For this reason, or for other 
reasons, we may expect serious difficulties in numerical computations 
based on (20) and (21) as they stand. These difficulties are considerably 
reduced by the transformations described in the next section. 

In the example under consideration there is of course no difficulty 
about inverting matrices analytically. It is found finally that 
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at epoch t 2 , 



where 



x = — __ {(r/[2(xi + -i ; 2) + t(vi + v 2 )] + <t v 2 t 2 x 2 \ , 

4l<t x ~ + T-a v - 



V = 



4(T X - + T-(X V 2 



and 



C = 



[2<r/(t;i + Vo + ar) + jr.Vfas - Xi + £ar 2 ) } , 



2<T- + T>, T(T 



4<r» 2 + rW 2 



to - ,, 



5<T« 



Under the assumptions implicit in Ci and C 2 , that xi , vi , :r 2 , and v 2 

are not correlated, the equations for x and v may be verified by a more 

familiar method. We have in fact two independent estimates of the 

position at epoch t 2 , viz., 

i - - 

Xi + — — j: — t with variance a x + — — , 

and 

x 2 with variance <r s . 

Taking the weighted average of these two estimates, each estimate being 
weighted in inverse proportion to its variance, we get x, as expressed 
above. The variance of x is the harmonic mean of the two variances, so 
that 

2/o 2 , 2 2\ 

var .i' = — j — — : — — j — , 
4a x - + r-a v - 

in agreement with the equation for C. 

Similarly, we have three independent estimates of the velocity at epoch 

t 2 , viz., 

V\ + (xt with variance a v , 
v 2 with variance o-„", 



and 



a,-o — Xi . 1 ... 2o- 2 

— + - a,T with variance — 7 

T 2 T- 
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Taking the weighted average of these three estimates, each estimate 
being weighted in inverse proportion to its variance, we get v, as expressed 
above. The variance of v is the harmonic mean of the three variances, so 

that 

var{5) = 2er.V/(W + iV), 

in agreement with the equation for C. 

The covariance of x and may also be verified, but the algebra is more 
involved. 

4.2 Transformation of Equal ions 

Omitting the explicit notational reference to updating (i.e., extrapola- 
tion) in (20) and (21), we have 

i= C.(C, '-ft + C.- 1 .*), (23) 

where 

C - (CV' + CV 1 )- 1 . (24) 

These equations require the inversion of three matrices which, at least 
in the case of angles-only data, are usually extremely ill-conditioned. 
Hence, they have been transformed in order to reduce this difficulty. 
The transformations described here were developed by A. J. Claus and 
the author. 
Introducing the identity 

Cr'-e, = id' 1 + Cf l )-t x - 6V'. ei 
into (23), we get 

e = e, - C-C 2 ~ I -(e, - fe). 
Now, 

c-cv 1 = (cr l + cr 1 )- 1 ^, -1 = [ft-ccr 1 + cr 1 )] -1 
= (i H-cvcrT 1 = [(ft + (?,). erf' 
= Ci-cci + cr 1 . 

Next, let 

/ J , = S-Ci-iS, P 2 = £-C 2 -£, (25) 

where S is a diagonal matrix in which each diagonal term is the reciprocal 
of the square root of the sum of the corresponding diagonal terms of 
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Ci and C 2 . (Thus, every diagonal term of Pi + P 2 is unity.) Then, 

Ci-(Ci + C 2 )-' = flT*-Pi.(Pi + Pj)" 1 -^ 
Hence, 

e = ei - flT'.Pj-fPi + Pir'-S-Cft - «). 

Similarly, 



= e 2 + S _1 -P 2 -(Pi + Pa) 1 -iS-(«i - €2). 



Next, let IFi and ll' 2 be arbitrary square matrices whose sum is a unity 
matrix. Then, 

- e = Wi-ei+ ll'.-eo - ft-(e, - e 2 ), 

where 

fi = (TTi-iSr'.p, - WrST l -Pt)'(Pi + p 2 ) -1 -s. 

Since 1\ + P 2 may yet he ill conditioned for inversion, we now use a 
well-known artifice, and write 

R = [WriT^iPx-G) - HV»S-'-(P,-tt)H(P, + PO-GT'-fl, 

where G' may he regarded as another arbitrary square matrix although 
it merely represents a set of rules for combining rows and/or columns of 
p x _j_ p 2t as we ll as of Pi and P 2 individually. [The normalization of the 
matrix sum (\ + C 2 to Pi + P 2 and the preservation of its symmetry 
by the introduction of the matrix S, as in (25), simplifies the implementa- 
tion of the matrix G as a set of operational rules.] Finally, Wi and W 2 
are restricted to diagonal matrices, so that 

R = S- l Wi'(PvO) - HV(P.-<?)]-[(Pi + P^-Gf'-S. (27) 

Noting that the right-hand member of (23) reduces to 2C if «i is re- 
placed by Ci and e 2 is replaced by C 2 , it follows from (26) that 

C - IPFi-Ci + WVCi - fl-(C, - d)]. (28) 

Tor further details, sec Rcf . 12. 

4.3 Debiasing Single-Pass Estimates 

Depending upon the method used to obtain single-pass estimates of 
the orbital elements in the first (or intrapass) stage of the pass-by-pass 
method of orbit refinement described in Section IV, the single-pass 
estimates may be biased on the average even if the errors in the observa- 
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tional data are not biased. Unless these biases are tolerable, they must 
of course be removed. In this section we will describe a method of remov- 
ing these biases in the single-pass estimates of the orbital elements. 
However, biases in the single-pass estimates of the orbital elements due 
to biased errors in the observational data will not be removed by this 
method. 

In the interest of simplicity, the description will be by analogy with 
a much simpler problem which involves only one observable coordinate, 
one parameter to be estimated, and does not involve time at all. We will 
consider two methods of estimation, of which the first is analogous to the 
classical method of orbit refinement, and the second is analogous to any 
method of obtaining single-pass estimates which are biased on the aver- 
age even if the errors in the observational data are not biased. 

Let y be a function of the observable coordinate x, and let 

x% = -To + a (i = 1, 2, ••• ,n) 

be the observed values of a:, where the e's are uncorrelated random er- 
rors with ave { e t j =0 and var { e,} = a 2 for every i. The most direct way 
of estimating y = ?/(.r ) is obviously to compute first 

x = {\/n)Y.Xi, 

and then y(x). The nature of the estimate obtained in this way is deter- 
mined as follows. Since 

•i"' = -l'u + «, 

where 

a = (1/71)2 «,-, 

then, to second-order terms in the e,'s, 

y(x) = y + a {l a + \b Q a, 
where a = dyo/dxo , and b — d'yu/dxo". Now, 

ave {a} =0, ave \a\ = a/n. 
Hence, to second-order terms in a, 

ave [y(x)\ = y„ + (boa/2n), 

and 

var {y(x)\ = ave {[y(x)f\ - [ave \y(x)\] 2 = a£a/n. 
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If a is not known, an estimate of var \y(x)\ is given by 

2 



Since 

lim ave {y(x)\ = y , 

n-»oo 

the estimate is asymptotically unbiased; and since 

lim var \y(x)\ = 0, 

H-»CO 

the estimate is "consistent" in the sense that the probability is unity 
that it will be correct to at least the second order as n — * <» . These re- 
sults are indicative of the nature of the estimates of orbital elements 
obtained by the classical method of orbit refinement, in which the esti- 
mates are such as to "predict" values of the observable coordinates 
which agree with the actual observations in the least squares sense. 

Now, consider another way of estimating y Q . We compute y(x) for 
each observed value of x, and define the estimate of y as 

?/„ = {\/n)YsVi where //,• = y(x { ). 

The purpose of estimating ;/ in this way is to permit the estimation of 
var \y \ without using a Q . Thus, 

where 

R< = Vi - yo ■ 

The importance of this is that the analog of a in the computation of or- 
bital elements from, say, two complete radar fixes (each fix consisting of 
a range and two angles) is the inverse of a 6 by 6 matrix whose compo- 
nents are functions of the epochs of the two fixes. The computation of 
this matrix may be avoided by computing a set of orbital elements from 
each pair of radar fixes, averaging over the sets, and estimating the co- 
variance matrix from the residuals. 

The nature of the estimate y is determined as follows. Since, to sec- 
ond-order terms, 

iji = 2/o + aoe. + ^''oe, , 
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then 













j/o = .'/.. + (ha + |&o/3, 




where a 


is 


as 


previously defined, and 














p = (i/»)£«. 




Now, 










ave j/3| = a'. 




Hence, 










ave ij/„j = |/o + (6 <r 2 /2), 




and 




















var 


f?7o! 


= ave {|(/ J J | - [ave |?/„|| 2 = 


2 2 / 

= an a fn. 



Thus, #„ is a biased estimate of //„ . In particular, it should be noted that, 
while the variance of this estimate decreases with increasing n, the bias 
in the estimate is independent of //. Hence, the variance bears no rela- 
tion to the accuracy of the estimate. 

The bias may be removed by the following supplementary procedure. 

1. Compute A o such that 

.y(.fo) = i/o . 

This is analogous to computing artificial tracking data (at the same 
epochs as the actual tracking data ) from the biased single-pass estimates 
of the orbital elements analogous to fju . The method of computing track- 
ing data from orbital elements must of course be numerically compatible 
with the method of computing orbital elements from tracking data in 
the absence of observational errors. 

2. Compute 

.r, = 2.fn — Xi , i = 1, 2, • •• ,n. 

This is analogous to combining the actual tracking data with the arti- 
ficial tracking data computed in the preceding step. (The choice of a 
combination such that the random error in each x ,■ is equal in magnitude 
but opposite in sign to the random error in the corresponding ;r, was 
suggested by D. R. Brillinger.) 

3. Compute 

Vi = y(&i), i = 1, 2, ••• , n. 
This is exactly the same procedure used in computing //, = y(Xi). 
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4. Compute 

Vi= *(3y«- in), i = 1,2, ••■ ,n. 

5. Compute 

£ = ( i/-/j ) J^Vi as f* 6 estimate of yo . 
The nature of the estimate y is determined as follows. To second-order 
terms, 

x = -to + a " J-*- (« 2 - /3) • 



Hence, 



6 = xo - (« ~ 2a) - - (« 2 - /3), t = 1, 2, • • •, n, 
2/,- = z/o - a„(«< - 2a) + £ (e* - 4a* + 2a« + 2/3), 



Iji = 2/o + Oo(2e,- - a) + ^ («< 2 + 2ae, - a - &), 



2/o = 2/o + aoa + i&oa • 
Thus, to second-order terms in the e/s, 2/0 is the same as y(x) . It is asymp- 
totically unbiased, and it is a consistent estimate of 2/0 . It may be abso- 
lutely debiased, to second-order terms in the e t % by changing step 4 of 
the supplementary procedure to compute 

Then, in determining the nature of the estimate 2/0 , we now have 
fH = ?/o + -S-r I(2n - 1) u - na] 

lb J- 

+ , b ° [(n - 1) u + 2nati - na - n$\, 
2(n — 1) 

27o = 2/o + ao« + 57 — — pj (na" - 0). 

Hence, 

ave {?7o| = 2/0 • 
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Generally, if step 2 of the supplementary procedure is changed to 
compute 

Xi = §[(1 + w)£o+ (1 - w)xi], 
and step 4 is changed to compute 

Vi = M(l + W) Vi + (1 - \V)y,], 
then, 

whence. 

Hence, the estimate ?y n is asymptotically unbiased if 

W = w * ± 7 
w 2 - 1 ' 

and is absolutely unbiased if 

2 , 7n + 1 
w + r 

ir = u " x 



w 2 — 1 

The choice of w should be made with some regard to the fact that 

(n - 1) (1 + w) (3 - w) 



var [Xi\ — a 



1 - 



An 



l 



The choice is iv = 3 in step 2 of the supplementary procedure, and 
W = 2 in step 4, so that var {£ ,| = a = var \x { \ . 

V. CLAUS'S METHOD 

The pass-by-pass method described in Section IV was used at the An- 
dover, Maine, station of the Telstar I experimental satellite communica- 
tions system. The intrapass estimates of the orbital elements were 
computed from angles-only data by a method described in some detail 
in Refs. 11 and 12. Suffice it here to say that a set of orbital elements is 
computed from each set of four sightlines (of which there may be as 
many as 200 in a single pass), the sets of orbital elements are averaged, 
and an estimate of the covariance matrix is computed from the residuals. 
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(Compare this outline with that of the example cited in Section 4.3, which 
uses complete radar fixes, including range. It should be noted that unless 
the number of sets of orbital elements is at least equal to the number of 
orbital elements the covariancc matrix will be singular.) This method 
gave excellent results for about four weeks. After that period it began to 
give sporadically bad results — typical errors of 10 5 feet in single-pass 
estimates of the semimajor axis. 

Extensive simulation studies by W. C. Ridgway III showed that most 
of the trouble was probably due to the sensitivity of single-pass estimates 
of the orbital elements to a bias error in the elevation angles, where these 
estimates are derived from single-tracker angles-only data. This result 
was subsequently confirmed by some formal analysis by the author. 
Ridgway's studies showed, in particular, that the sensitivity increases 
rapidly with decreasing length of pass, and with increasing maximum 
elevation angle, although sightlines at elevation angles over 82.5 degrees 
(or under 7.5 degrees) were not used. This is consistent with the fact 
that the sporadically bad results began to occur when the perigee of 
Telstar I had precessed sufficiently to make the passes at Andover sub- 
stantially shorter than they were during the first week, and a substan- 
tially larger proportion of the passes had high maximum elevation angles. 
Ridgway's studies showed that, under these conditions, an error of 10 5 
feet in the single-pass estimate of the semimajor axis could easily be due 
to a bias of 0.01 degree in the elevation angles. 

In order to overcome the sensitivity of a single-pass single-tracker 
angles-only method to a bias error in the elevation angles, on short passes 
with high maximum elevation angles, A. J. Claus has developed a method 
of orbit refinement which permits the use of a few, perhaps only one or 
two, measurements of range in each pass. This method was intended 
to be used in the intrapass stage of the pass-by-pass method described in 
Section IV, but it may be used as a self-sufficient method, just as Swerl- 
ing's method may be used either in the intrapass stage of the pass-by- 
pass method or as a self-sufficient method. 

Although Claus developed his method with no foreknowledge of 
Swerling's method, it turns out that his method differs from Swerling's 
method only in the explicit introduction of an iterative routine which, 
as in the classical method, improves the estimates of the orbital elements. 

Instead of (12) we now substitute (1) into (11) so that 

Q= (e-^'-tfrMe-*) (29) 

+ [r - J-(e - e„)]'-3> -1 -[r - J-(e - «.)], 
where J and r are defined by (2) and (3), and e has the same signifi- 
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cance as in the classical method. Now, Q is minimum with respect to e 
if 

Ci~ l -(* - «i) - J'-* _1 -[r - ./-(e - eo)] = 0. 
Denoting the value of e which satisfies this equation by e, we have 

I = C'[C l " 1 -h+ ./'^"'-./•60+p], (30) 

where 

C = [&->+ ./'.*-'. J]~\ (31) 

p = J'.qT l r. (32) 

Equation (30) may be expressed in the form 

e = ei + ('•[./'•* ■'../•( e„ - I,) + p], (33) 

where C, ./, and p (the latter through r as well as ./) depend implicitly 
on e,i . Comparing this with (15), it will be seen that the iterative use of 
(30) or (33), substituting l for e at each iteration, until the difference 
between e and I is negligible, is precluded in Swerling's method by the 
constraint e = ei . 

Finally, since (29) may be expressed in the form 

Q = (e — e)'-C"'-(e — i) + terms independent of e, (34) 

it follows that C is the covariance matrix of e. 

vi. battin's method 

This method, as described in Ref. 14, is essentially a special case of 
Swerling's method. However, by introducing new data one at a time, 
Battin's method avoids the inversion of matrices. (Swerling's JAS article 
contains equations whereby matrix inversions are avoided if new data 
are introduced one at a time, but these equations appear at the end of the 
section entitled "Statistics of Propagated Errors," and their use for 
avoiding matrix inversions is not explicitly stated.) 

In the notation of Section III, (15) may be written in the form 

6=^+ [J'-.J + SCrV-J'-r, (35) 

where a is the variance of the scalar $, ./ is a 1 by matrix (i.e., ./' is a 
6-rowed vector), and r is a scalar. Now, if 

a = ./■(',■./' + a 1 (a scalar), (36) 

then, 
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_, ., JC yJ' + cr 2 
./ — J " - > 



= -[f-J&+ A- J' , 

a 



(37) 



Substituting (37) for the last /' in (35), we get 

1 -* + £&- j\ (38) 

a 

Further, since (16) may be written in the form 

C = a%J'-J + oVJfT" 1 , 

then, 

c- & = aj'-j + '&'T 1 - a, 

and, therefore, 

[j'.J + «- f ft-*]-(C- ft) = -J'JCi. 
Substituting (37) for ./' in the right-hand member, we get 

C = Cr --CrJ'J-Ci- (39) 

a 

Equations (36), (38), and (39) are equivalent to equations (30), (29), 
and (33) in Battin's paper (Ref. 14). 

From this description of Battin's method, it is evident that the in- 
version of matrices may be avoided also in Claus's method if new data 
(perhaps only the range data) are introduced one at a time. Comparing 
(33), (31), and (32) with (15), (16), and (17), it is clear that (38) and 
(39), with (36), are valid in Claus's method if 

r = j>-<p(e ) + J- (6o- «i), (40) 

and 

J = dtpM/deo. (41) 

The initial value of e may be taken equal to h , but thereafter i is re- 
peatedly substituted for e until the difference between I and « is negli- 
gible. 
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VII. MOD-S AND MOD-C METHODS 

Swerling's method and Claus's method can be modified to permit the. 
introduction of new data n at a time, where n < 6, at the cost of having 
to invert a matrix of order n. This might be worthwhile for n = 2 or 3. 

With regard to Swerling's method, let 

A = /•<?!•/'+* (an n by n matrix) . (42) 



Then, 



J' .<*>-' = J'.$- l .[J.C x -J' + ^)-A~ l 
= [/'•#"* J + 6r x ]-C v J'A 



(43) 



Substituting this into (17), substituting the resultant expression for p 
into (15), and taking account of (10), we get 

e = e, + CrJ'-.r'-r. (44) 

Equation (36) is a special case of (42), (38) is a special case of (44), and 
(39) is a special case of 

C = C, - CvJ'A-'-J-d. (45) 

With regard to Claus's method, (42), (44), and (45) are valid if /• 
and ./ are defined by (40) and (41 ), and the repeated substitution of 
e for to is carried out until the difference between c and i is negligible. 

VIII. SUMMARY AND CLOSING REMARKS 

8.1 Summary 

The classical "differential corrections" method of orbit refinement, 
occasionally supplemented by the use of "normal places," is appropriate 
for astronomical bodies whose relative positions change very slowly, 
whose relative angular positions, viewed from the earth, can be measured 
with extreme optical precision, and which therefore require compara- 
tively small quantities of observational data to establish their orbits with 
great accuracy. However, that method is very unwieldy for artificial 
earth satellites or short-range space probes, where the relative inaccuracy 
of the observational data must be offset by greater quantities of observa- 
tional data. 

For artificial earth satellites or short-range space probes, Swerling's 
method is more practical, chiefly because it does not require all of the 
observational data to be processed together. However, to omit any part 
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of the old observational data which has been processed by Swerling's 
method it is necessary to start all over again in the sense that all of the 
old observational data which are to be retained must be reprocessed in 
the same way, along with any new observational data which might be 

available. 

The pass-by-pass method, on the other hand, operates on the princi- 
ple of computing an independent set of estimates of the orbital elements 
for each pass, and combining the sets of single-pass estimates in an opti- 
mum way. Thus, entire blocks of observational data may be omitted 
without actually reprocessing any of the old observational data from 
which single-pass estimates of the orbital elements have already been 
computed. The method of obtaining the single-pass (intrapass) estimates 
is optional. In the Telstar I experiments, it consisted in dividing up the 
data (no range data) into mutually interlaced independent sets of four 
sightlines, computing a set of orbital elements from each set of four 
sightlines, averaging over the sets of orbital elements, and computing a 
covariance matrix for the average set, from the residuals. This particular 
intrapass method introduces biases (apart from the biases in the data) 
and a method of eliminating these computational biases was developed 
but was not used in the Telstar I experiments. After four weeks of suc- 
cessful operation, a more serious source of trouble arose, which was 
traced to the increasing sensitivity to bias (residual boresight error and 
sample bias) in the elevation angle data. This increasing sensitivity was 
associated with the precession of perigee to latitudes close to that of the 
tracker. 

To overcome the sensitivity of the angles-only intrapass method to 
bias in the elevation angle data, Claus developed a method, intended to 
be used chiefly in the intrapass stage of the pass-by-pass method for 
Telstar II, which could accept occasional range data. This method is 
essentially Swerling's method supplemented by an iterative routine 
which improves its accuracy. 

In Swerling's or Claus's method, the inversion of six-by-six or higher- 
order matrices can be avoided by borrowing a detail from Battin's 
method of spacecraft navigation, provided that the observational data 
are processed only one at a time, as in Battin's method. However, at the 
cost of inverting n-by-n matrices, where n < 6 (say, 2 or 3) , the observa- 
tional data may be processed n at a time, by modified forms of Swerling's 
or Claus's method. 

The comparisons of the methods described in this paper have been 
made only to the extent that motivated the development of the newer 
methods. The practical details of these methods are interchangeable to 
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a large extent, so that different and more appropriate combinations of 
these details may be made for other specific practical applications. 

8.2 Closing Remarks 

Clans has pointed out that, in principle, none of the alternative meth- 
ods described in this paper, including the use of normal places in the 
classical method, can be as efficient, in a statistical average sense, as the 
classical method without normal places. The reason for this is simply 
that the classical method without normal places allows the maximum 
possible freedom in fitting the estimates of the orbital elements to the 
observational data. Hence, the choice of a method, or combination of 
methods, for a particular application, usually involves a small sacrifice 
in accuracy. 
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APPENDIX 

Equations (20) and (21) may be derived in another way (essentially 
that used by Aitken in Ref. 3) which does not depend upon the minimiza- 
tion of a quadratic form and provides further insight into the signifi- 
cance of these equations. 

Let x be the true value of an >?-rowed vector (one-column matrix). 
Let .f be an unbiased estimate of .r, with 

avej(.r - x)-(x - x)'\ = C, 

where the prime stands for transposition and "ave" stands for ensemble 
average. The nth-order square matrix C is the covariance matrix of .? . 
Let x be another unbiased estimate of .r, with 

ave{(.f - x)-(x - x)'\ - C, 

and let it be assumed that .r and x are independent, so that 

ave {(.r - x)-(x - .r)'! = 0. 

Now, consider the weighted linear average 

x = ir-.T + W-x (46) 
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where W and W are nth-order square matrices, with 

17 + W = I (47) 

where / is the nth-order unit matrix. Since 

x - x = W-(x - x) + W-(x - x), 
it readily follows that if 

C = ave [(z — x) -(x — x)'\ 
then 

C = WCW' + W-CW'. (48) 

Each of the diagonal elements of C, viz., 

Cu = Z (ftii'Gjk-ffi* + Wij-Cjk-Wnc) 

i.k 

will be minimized under constraints on W ik + W« by minimizing 
<?«- 2 £\« ■(#,•* + #<*) 

k 

where the \«'s are Lagrange multipliers. This requires that 

3 J 

for every i, k. In matrix notation, 

W-C = X, JF-C = X 
where X is the nth-order square matrix of the X,*'s. Hence, 

W = \-C~\ W = \-C~\ (49) 

where, to satisfy (47), 

X = (C'- 1 + C- 1 )-. (50) 

By (48) and (49) 

C = X-(C ,_1 + &*)+', 
= X'by (50), 
but, since X is a symmetrical matrix, 

C = X. (51) 

Finally, by (46), (49), and (51), 

x = C-(C'- l -x + CT'-x), (52) 
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where, by (50) and (51), 

C = (C- 1 + C-)- 1 . (53) 

This derivation shows that the diagonal elements of C, which are the 
variances of the components of .f;, have been minimized independently 
of one another. 
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